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ABSTRACT 


Paacal-SC  (Pascal  for  Scientific  Computation)  is  an  extension  of  the 
popular  computer  language  Pascal.  Pascal-SC  greatly  facilitates  scientific, 
engineering,  and  statistical  calculations  on  microcomputers.  The  most  important 
of  the  additional  capabilities  of  Pascal-SC  arel  v~ 

\ 

'—ttf-  Accurate  floating-point  arithmetic  for  real,  complex,  and  interval 
numbers,  vectors,  and  matrices/-"  Expressions  involving  these  types  can  be 
written  using  essentially  ordinary  mathematical  notation,  making  programming  and 
documentation  easier.  Numerical  results  are  generally  obtained  with  more 
accuracy  than  with  conventional  floating-point  arithmetic;  in  particular,  scalar 
products  of  vectors  are  competed  to  the  closest  floating-point  number.  Standard 
routines  are  provided  which  return  guaranteed  error  bounds  as  well  as  answers 
for  the  solution  of  linear  systems  of  equations,  matrix  inversion,  and 
eigenvalue-vector  calculations. 

c.XTZ"'~ . 

~  i-Wrl  User-defined  operators  to  allow  the  manipulation  of  data  of  various 
nonstandard  types  by  expressions  written  in  essentially  ordinary  mathematical 
notation,  instead  of  the  sequences  of  calls  to  functions  and  procedures  required 
in  ordinary  Pascal.,  This  capability  permits  simple  and  straightforward  use  in 
programs  of  various  coordinate  systems  for  representation  of  variables, 
arithmetic  for  polynomials  and  series,  automatic  differentiation,  etc.,  and  thus 
adds  considerable  flexibility  and  power  to  the  language. 

Brief  descriptions  of  these  useful  features  of  Pascal-SC  are  given, 
together  with  illustrative  examples.  Some  familiarity  with  ordinary  Pascal  is 
assumed.  . 
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SIGNIFICANCE  AND  EXPLANATION 


The  advent  of  the  personal  microcomputer  has  added  a  new  dimension  to 
scientific  computation.  Numerous  engineers,  scientists,  and  statisticians  are 
using  such  machines  at  work  and  at  home  to  solve  both  routine  and  research 
problems  which  arise  in  their  professions.  In  addition  to  providing  answers  to 
specific  problems,  microcomputers  can  be  used  to  try  various  methods  on  scaled- 
down  problems  in  order  to  determine  which  are  the  most  promising  for 
implementation  on  mainframe  systems.  A  crucial  factor  for  the  efficient  use  of 
microcomputers  is  the  programming  language  used.  This  report  describes  the 
language  Pascal-SC  (Pascal  for  Scientific  Computation),  which  offers  unique 
advantages  to  microcomputer  users. 

First,  Pascal-SC  employs  really  accurate  floating-point  arithmetic  for 
real,  complex,  and  interval  numbers,  vectors,  and  matrices.  Ordinary  operator 
notation,  rather  than  calls  to  functions  and  procedures,  makes  these  arithmetics 
easy  to  use.  Furthermore,  guaranteed  error  bounds  are  returned  by  the  library 
programs  for  solution  of  linear  systems  of  equations,  matrix  inversion,  and 
eigenvalue-eigenvector  calculations. 

Second,  the  user  can  define  operators  for  arbitrary  data  types  used  in  the 
program.  This  convenience  of  operator  notation  makes  programs  easier  to  write, 
understand,  and  document.  Furthermore,  the  resulting  programs  are  generally 
much  shorter  than  the  corresponding  programs  written  in  ordinary  Pascal. 

Finally,  Pascal-SC  is  a  simple  extension  of  Pascal,  a  widely  used  language 
for  microcomputers.  The  Pascal-SC  compiler  will  translate  programs  written  in 
ordinary  Pascal,  and  programmers  who  already  know  Pascal  can  learn  the  extra 
features  of  Pascal-SC  programming  in  a  few  minutes.  Thus,  none  of  the 
investment  in  learning  Pascal  or  in  programs  written  in  Pascal  is  lost  in  going 
to  Pascal-SC. 

This  report  explains  the  basic  features  of  Pascal-SC,  assuming  some 
familiarity  with  ordinary  Pascal,  and  illustrates  them  by  examples .  Although 
this  document  does  not  proport  to  be  a  user's  manual,  it  could  serve  as  a  guide 
to  a  programmer  already  proficient  in  Pascal.  Several  years  of  experience  with 
Pascal-SC  on  various  microcomputers,  including  the  Zilog  MCZ-1  at  the 
Mathematics  Research  Center,  verifies  that  this  language  is  a  flexible  and 
powerful  tool  for  research  into  numerical  methods,  as  well  as  for  routine 
calculations. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC ,  and  not  with  the  author  of  this  report. 


AN  INTRODUCTION  TO  THE  SCIENTIFIC  COMPUTING 


LANGUAGE  PASCAL -SC 


L.  B.  Rail 


1.  Objective.  A  brief  description  of  the  language  Pascal-SC  (Pascal  for  Scientific 
Computation)  will  be  given  to  explain  features  of  this  extension  of  ordinary  Pascal  (10) 
which  are  particularly  useful  for  scientific,  engineering,  and  statistical  calculations, 
particularly  on  microcomputers •  The  reader  is  assumed  to  be  familiar  with  Pascal 
programming  on  at  least  an  introductory  level,  and  to  have  had  some  experience  with 
numerical  computation  of  the  kind  which  arises  in  scientific  and  engineering  problems. 

With  this  background,  it  should  be  easy  to  appreciate  the  advantages  of  the  additional 
features  of  Pascal-SC. 

In  order  to  keep  the  discussion  short  and  to  the  point,  many  details  will  be  omitted, 
and  a  formal  description  of  the  language  will  not  be  given.  For  operational  and 
programming  details,  the  reader  should  consult  (3),  [19],  [28];  formalities  are  given  in 
[S].  Here,  simple  examples  will  be  used  to  illustrate  ideas  as  they  are  introduced. 

2.  Why  Pascal-SC?  As  pointed  out  by  Wirth  [10],  the  introduction  of  a  new  computer 
language  requires  careful  justification.  The  same  applies  to  an  extension  or  modification 
of  an  existing  language,  particularly  a  language  which  is  as  successful  and  widely  used  as 
Pascal.  The  most  important  additional  features  of  Pascal-SC  are: 

(i)  Accurate  floating-point  arithmetic  with  controllable  rounding; 

(11)  User-defined  operators  to  facilitate  progressing  and  documentation. 

Furthermore,  it  is  of  considerable  importance  to  note  that: 

(ill)  Pascal-SC  retains  the  features  of  ordinary  Pascal. 

i 
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Thus,  non*  of  th*  investment  in  learning  to  program  in  Pascal  or  in  programs  already 
written  in  Pascal  is  lost  in  going  from  Pascal  to  Paacal-SC.  The  Pascal-SC  compiler  can 
translate  programs  written  in  ordinary  Pascal.  Moreover,  programming  in  Pa seal -EC  will 
corns  vary  naturally  to  the  Pascal  programmer,  as  will  be  seen  from  the  examples  given 
below. 

With  regard  to  (i),  the  floating-point  arithmetic  provided  by  Pascal-SC  is  implemented 
not  only  for  real  floating-point  numbers  (type  REAL),  but  also  for  complex  numbers, 
intervals,  and  vector*  and  matrices  of  these  types.  This  allows  convenient  and  accurate 
computation  with  the  kinds  of  numerical  data  most  frequently  encountered  in  scientific  and 
engineering  problems.  The  floating-point  arithmetic  of  Pascal-SC  [28]  is  constructed  on 
the  basis  of  the  theory  of  computer  arithmetic  given  by  Kulisch  and  Miranker  [14],  which 
guarantees  accuracy,  controllability  and  reliability  of  th*  results.  In  order  to  keep  the 
coeg>iler  small  enough  to  be  convenient  to  use  on  microcomputers  and  still  provide  these 
additional  features,  extensive  use  is  made  of  external  libraries  of  pretranslated  code 
which  the  compiler  can  link  to  the  user's  program,  or  source  code  which  can  be  compiled  as 
part  of  it. 

The  second  important  additional  capability  of  Pascal-SC  is  that  it  allows  user-defined 
operators  to  permit  the  manipulation  of  nonstandard  data  types  in  ordinary  mathematical 
notation.  For  example,  if  A  is  a  matrix,  end  x,  b,  c  are  vectors,  then  the  programmer  cen 
write  the  statement 

(2.1)  c  A*x  +  bi 

in  Paacal-SC  to  perform  the  indicated  calculations.  This  notation  follows  ordinary 
mathematical  usage,  and  thus  has  th*  advantages  of  clarity  and  simplicity  compared  to  the 
calling  of  procedures  and  functions  to  obtain  th*  same  result  in  ordinary  Pascal.  In  order 
for  th*  Peecal-SC  compiler  to  accept  (2.1),  the  heading  of  th*  program  has  to  contain  a 
definition  of  the  binary  operators  *  to  perform  matrix  by  vector  multiplication  and  +  to 
perform  vector  addition.  (Source  code  for  these  operetora  is  included  in  th*  Paseal-SC 


package  for  vector  and  matrix  arithmetic.)  In  addition  to  this  "overloading*  of  standard 
operator  symbols,  Pascal-SC  permits  the  user  to  give  operators  arbitrary  names  (for 
example,  XOR  for  "exclusive  or"),  and  assign  priorities  to  such  operators.  One  particular 
convenience  of  Pascal-SC  is  that  the  operator  **  can  be  defined  to  perform  exponentiation, 
which  makes  the  writing  of  polynomials  and  other  functions  containing  powers  simpler  than 
in  ordinary  Pascal.  In  allowing  usar-dafined  operators,  Pascal-SC  is  similar  to  Algol  68 
and  Ada. 

These  points  will  be  discussed  in  more  detail  in  the  following  sections. 

3.  Floating-point  REAL  arithmetic.  This  is  the  "built-in*  arithmetic  of  Pascal-SC 
for  floating-point  numbers  (type  RZAL).  Since  this  arithmetic  is  based  on  the  general 
theory  of  computer  arithmetic  given  in  114],  it  is  accurate,  controllable,  and  reliable. 
Before  going  on  to  details,  a  precise  statement  of  the  meaning  of  these  terms  is  necessary 
because  the  related  concepts  of  "accuracy"  (the  exactness  with  which  results  are 
calculated)  and  "precision"  (the  number  of  digits  used  in  the  representation  of  floating¬ 
point  numbers)  are  often  confused.  For  example,  the  result  32.0  -  31.0  »  1.00  is 
calculated  with  low  precision  (3  decimal  digits),  but  high  accuracy  (exactly).  By 
contrast,  unvac  1100  floating-point  hardware  gives 


134217728.0  -  134217727.0  «  2.00000000, 


which  is  done  with  higher  precision  (9  significant  digits),  but  no  accuracy,  since  the 
correct  answer  is  1.00000000  [21]. 

It  is  possible  to  discuss  the  idea  of  accuracy  Independently  of  the  particular 
precision  used,  since  each  floating-point  system  contains  only  a  finite  set  of  numbers.  In 
a  given  system  8,  two  floating-point  numbers  u,v  with  u  <  v  will  be  said  to  be  adjacent  if 
there  is  no  floating-point  number  w  such  that  u  <  w  <  v.  For  x,y  e  S,  the  exact  result  xOy 
of  an  arithmetic  operation  o,  where  o  €  {♦,  -,  *,  /},  will  either  be  a  floating-point 
number,  or  a  real  number  w  such  that  u  <  w  <  v,  where  u,v  are  adjacent  floating-point 
numbers.  In  order  to  produce  a  floating-point  number  in  the  latter  case,  w  is  "rounded"  to 


•n  element  of  S,  which  should  bs  either  u  or  v.  This  is  the  basic  requirement  for 
reliability  of  a  floating-point  arithmetic  operation.  The  floating-point  software  supplied 
for  some  microcomputers  does  not  mest  this  simple  requirement,  and  neither  does  the 
floating-point  hardware  of  some  mainframe  computers,  nils  unhealthy  situation  has  lead  the 
IBSE  to  undertake  the  promulgation  of  standards  for  floating-point  arithmetic  (see  the 
SIGNUM  Newsletter  for  October,  1979). 

An  accurate  rounding  R  of  the  floatingpoint  operation  o  selects  R(xPy)  equal  to  u  or 
v  to  minimite  the  roundoff  error 

(3.1)  e(x<5y)  -  | R(xOy)  -  (xOy )  |  . 

This  is  the  beet  possible  answer  (bfa)  rounding  [29].  In  case  of  a  tie  between  u  and  v,  a 
suitable  rule  is  invoked.  In  Pascal-6C,  this  rounding  is  to  the  one  of  u,v  which  is 
furthest  from  sero,  in  order  to  satisfy  the  condition  of  antisymmetry  of  R,  R(-x)  -  -  R(x), 
which  is  required  by  the  general  theory  [14]. 

The  distance  |u  -  v|  between  u  and  v  are  will  depend  on  the  precision  of  the  floating¬ 
point  numbers  being  used)  this  determines  the  maximum  roundoff  error  of  a  reliable 
calculation.  Of  course,  if  the  arithmetic  of  the  machine  being  used  is  not  reliable,  then 
roundoff  error  is  not  related  in  a  simple  way  to  precision,  and  the  attempt  to  "buy*  more 
accuracy  by  using  increased  precision  can  be  futile,  as  well  as  expensive. 

niere  are  other  ways  in  which  rounding  to  a  reliable  result  can  be  carried  out, 
including! 

(1)  x  is  rounded  upward  to  vi 
(ii)  x  is  rounded  downward  to  u. 

The  direotatf  roundings  (i)  and  (ii)  are  necessary  to  support  interval  arithmetic  [17], 
[16],  [14],  among  other  things.  Rounding  in  a  floating-point  arithmetic  is  said  to  be 
controlled  if  the  user  can  choose  the  method  desired  for  the  result  of  a  given  floating¬ 
point  arithmetic  operation.  Pascal-SC  gives  tha  user  the  choice  of  BRA  and  the  upward  and 
downward  directed  roundings,  which  maans  that  a  total  of  twelve  operators  are  provided  for 


-4- 


the  four  basic  arithmetic  operations  +,  -,  *,  /,  and  the  three  roundings  listed  above: 


+  -  *  / 

{  BPA  rounding  } 

+>  ->  *>  /> 

{  Upward  rounding  } 

+<  -<  +<  /< 

{  Downward  rounding  } 

Thus,  the  Pascal-SC  programmer  can  control  the  direction  of  rounding  if  desired,  for 
example,  to  obtain  guaranteed  lover  or  upper  bounds  for  the  values  of  arithmetic 
expressiona  [6].  It  also  follows  from  the  reliability  of  Pascal-SC  arithmetic  that  the 
addition  and  multiplication  operators  (3,2)  are  commutative,  which  is  not  true  for  the  kind 
of  floating-point  arithmetic  ordinarily  encountered. 

For  the  microcomputer  implementations  of  Pascal-SC,  decimal  arithmetic  is  used,  and 
the  precision  of  floating-point  numbers  is  twelve  decimal  digits  in  scientific  notation, 
with  an  exponent  range  in  powers  of  10  from  -99  to  +99.  The  smallest  and  largest  positive 
numbers  are  thus  (UNREAL  -  1.00000000000E-99  (-  10_")  and  MAXREAL  -  9.99999999999E+99, 
respectively.  Zero  is  represented  by  0  *  O.OOOOOOOOOOOE+OO,  as  usual  [28].  The  use  of 
decimal  arithmetic  avoids  the  errors  introduced  by  conversion  between  binary  and  decimal 
upon  input  and  output.  Decimal  values  are  represented  internally  by  two  BCD  digits  per 
byte.  The  precision  of  the  Pascal-SC  floating-point  number  system  for  microcomputers  is 
thus  adequate  for  the  representation  of  most  numerical  quantities  of  interest  in  scientific 
and  engineering  computation.  Furthermore,  since  Pascal-SC  floating-point  arithmetic  is 
reliable,  accurate,  and  controllable,  there  is  seldom  any  need  for  more  than  12  decimal 
digits  of  precision  in  a  given  computation. 

To  illustrate  the  features  of  Pascal-SC  arithmetic,  consider  the  product 

(3.3)  A  -  (13.4565432278X0.000453782392145), 

The  emmet  value  of  A  is  not  a  12-digit  floating-point  number,  so  rounding  will  take  place 
in  the  corresponding  floating-point  multiplication  operations.  The  floating-point 
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operation  *  gives 


(3.4)  B  -  1 . 34565432278E+0 1  •  4. 53782392145E-04  -  6.  1063423759 IE-03, 

as  the  BPA  for  A,  which  is  in  error  by  at  most  -  5.0x10“15,  since  Paecal-SC  arithmetic  is 
reliable  and  the  distance  between  B  and  its  two  neighboring  floating-point  numbers  is 
4  ■  10  More  precisely,  (3.4)  establishes  that  A  belongs  to  the  half-open  interval 

(B  -  4,  B  +  j) ,  since  ties  are  rounded  away  from  0.  The  operations  *<  and  •>  with  directed 
rounding  give 


C  -  1.34565432278E+01  *<  4. S3782392145E-04  -  6. 10634237591E-03, 

(3.5) 

D  -  1 . 345654322788+0 1  *>  4. 53782392145E-04  -  6. 10634237592E-03, 

respectively.  Since  C  «  B,  the  result  (3.5)  shows  that  the  exact  answer  A  belongs  to  the 
interval  IC,D ]  “  lB,D]  <*  [8,  8  *  5);  therefore,  on  the  basis  of  (3.4),  a  belongs  to  the 
half-open  interval  (B  -  B  +  |)  n  (B,  B  +  5]  -  (B,  B  +  -|)  (see  Figure  3.1).  Thus,  this 
calculation  proves  that  A  j>  B  and  A  -  B  <  5. 0* 10-15.  This  gives  a  more  accurate  location 
for  A  than  the  BPA  answer  (3.4). 

A 

4 - E - ) - i - - 

B  -  -y  B-C  B+T  D  D+7 

Figure  3.1.  A  Detail  of  the  REAL  Floating-Point  Screen. 

In  addition  to  the  twelve  rounded  arithmetic  operations  listed  above,  Fascal-SC 
provides  the  most  frequently  used  standard  functions,  computed  to  BPA  accuracy  [28].  In 
order  to  keep  the  microcomputer  version  of  the  Pascal-SC  system  small,  extensive  use  is 
made  of  external  libraries,  so  the  compiler  will  bring  in  code  only  for  functions. 


procedures,  and  operators  actually  needed  by  the  program. 

5,  Floating-point  COMPLEX  arithmetic.  In  Pascal-SC,  manipulation  of  complex  numbers 
is  accomplished  by  subroutines  for  operators,  functions,  and  procedures  which  are  stored  in 
external  libraries.  For  this  reason,  the  declaration  of  complex  numbers  has  the 
stereotyped  form 

TYPE  COMPLEX  «  RECORD  RE , IM :  REAL  END; 

Thus,  the  representation  of  a  complex  number  Z  in  Pascal-SC  is  in  Cartesian  coordinates. 

For  input  and  output  of  Z,  the  standard  format  is  (Z.RE,  z.IM). 

In  ordinary  Pascal,  addition  and  other  arithmetic  operations  with  complex  numbers  have 
to  be  done  by  procedures  and  functions.  In  Pascal-SC,  however,  operator  overloading 
simplifies  notation  in  the  program  considerably.  To  Illustrate  this,  consider  as  a  simple 
example  the  source  code  to  enable  the  operator  +  to  add  complex  numbers: 

OPERATOR  +  ( A,B:  COMPLEX)  RES:  COMPLEX; 

VAR  (7:  COMPLEX; 

BEGIN 

U.RE  :»  A. RE  +  B.RE; 
tJ.IM  :«  A.  IM  +  B.  IM; 

RES  :«  U 
END; 

The  actual  coding  for  this  operator  is  essentially  the  same  as  for  a  function  or 
procedure  for  the  same  purpose.  However,  to  add  the  two  complex  numbers  V,W  and  assign  the 
result  to  the  complex  number  Z,  one  writes  only 

(5.1)  Z  :-  V  +  W; 

in  the  subsequent  program.  Addition  operators  have  to  be  defined  for  all  pairs  of  operands 
of  types  INTEGER,  REAL,  and  COMPLEX,  since  these  would  occur  naturally  in  expressions  being 
evaluated.  If  X,  R,  C  denote  generic  variables  of  types  INTEGER,  REAL,  and  COMPLEX, 
respectively,  then  six  addition  operators  are  needed: 
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(5.2) 


+C,  It  ♦  C,  C  ♦  X,  R  +  C,  C  +  R,  C  +  C. 


Similarly,  six  subtraction  operators 


(5.3)  -C,  X  -  C,  C  -  X,  R  -  C,  C  -  R,  C  -  C, 


are  required,  aa  well  as  five  multiplication  and  five  division  operators; 


(5.4) 


K*C,  C*K,  R*C,  C*R,  C*C, 


K/C,  C/K,  R/C,  C/R,  C/C. 


All  22  of  the  operators  { 5.2)— ( 5.4)  are  provided  in  an  external  library  in  the  form  of 
pretranslated  code,  and  the  corresponding  declarations,  for  example. 


OPERATOR  +  (RA:  REAL/  ft  COMPLEX )  RES:  COMPLEX; 
EXTERNAL  155; 


are  available  to  the  programmer  in  an  external  text  file.  Actual  coding  is  also  simplified 
considerably  the  fact  that  the  programmer  can  direct  the  compiler  to  refer  to  external 
libraries  for  type  declarations  and  definitions  of  operators  and  other  needed  functions  and 
procedures  (19),  (28).  The  use  of  this  feature  makes  the  source  code  for  a  Pascal-SC 
program  more  compact  and  readable.  Examples  of  programs  using  such  directives  are  given  in 
Appendices  A  and  B. 

All  complex  floating-point  operations  in  Pascal-SC  calculate  the  real  and  imaginary 
parts  of  the  result  to  BPA  accuracy.  Por  addition  and  subtraction,  operators  similar  to 
the  one  given  above  for  complex  addition  are  satisfactory;  however,  multiplication  and 
division  require  special  algorithms  to  attain  this  accuracy  (14),  (28).  For  example, 


-8- 


consider  the  function  CPIV  which  does  complex  division  by  the  usual  formula: 


(5.5)  FUNCTION  CD IV ( A , B  :  COMPLEX)  :  COMPLEX; 

VAN  DEN DM :  REAL; 

U:  COMPLEX; 

BEGIN 

DENOM  :“  B.RE  *  B.RE  +  B.IM  *  B.IM; 

U.RE  s»  (A. RE  *  B.RE  +  A.IM  *  B.IM)  /  DENOM; 
U.IM  :«  (A.IM  *  B.RE  -  A. RE  *  B.IM)  /  DENOM; 
CDIV  :-  U 
END; 


The  Pascal~SC  operator  for  complex  division, 


(5.6)  OPERATOR  /  (A,B:  COMPLEX)  RES:  COMPLEa; 

EXTERNAL  182; 


has  a  number  of  advantages  over  the  function  CPIV.  Among  these  are: 

a.  Accuracy.  Hie  results  of  Pascal-SC  complex  operations  are  calculated  to  BPA 
accuracy  in  the  sense  that  their  real  and  imaginary  parts  are  given  to  BPA  accuracy. 
Ordinary  Pascal  functions  and  procedures  for  COMPLEX  multiplication  or  division,  such  as 
CDIV  given  above,  cannot  attain  this  accuracy  because  of  the  number  of  roundoff  errors 
which  occur.  In  fact,  catastrophic  cancellations  can  happen  which  render  the  results 
almost  meaningless  when  the  ordinary  formulas  in  (5.5)  are  used.  For  example,  for 


V-  (1.23456789,  1.23456789), 


(5.9) 


W  -  (1. 0000123 E“ 05,  1.0000321E-05), 


-9- 


on*  9*t* 


(5.10) 


V/W  -  ( 1 . 23454048308E+05,  -1 . 222 16794612E+00 ) , 


CDIV(V,W)  -  (1.23454048308E+05,  -1.22216773503E+00 > , 


where  the  incorrect  digit*  of  CDIV(U,V)  are  indicated  in  boldface.  Even  in  thia  fairly 
harmless-looking  case,  the  algorithm  (5.5)  has  lost  five  significant  digits  in  the 
imaginary  part  of  the  quotient  in  a  single  division,  while  all  the  digits  of  the  Pas cal -SC 
result  for  V/W  are  correct.  To  be  sure,  roundoff  error  will  also  increase  with  repeated 
use  of  the  Pascal-SC  division  operator  (5.6),  but  at  a  alower  and  more  predictable  rate. 

b.  Deferred  overflow.  In  the  Paacal-SC  algorithms  for  complex  multiplication  and 
division,  overflow  does  not  result  unless  the  real  or  imaginary  part  of  the  result  is  > 

MAX REAL  in  absolute  value,  whereas  overflow  can  occur  in  (5.6)  in  the  calculation  of  ths 
Intermediate  values  DEN OH,  U.RE,  0.1M,  even  though  the  actual  result  has  real  and  imaginary 
parts  which  are  representable  by  floating-point  numbers.  For  example,  for 


(5.11)  V  -  (3.0E+99,  -1.0E+99),  V  -  (1.0E+99,  -1.0E+99), 


Pascal-SC  complex  division  gives  the  result 


(5.12)  V/W  -  ( 2 . 00000 00000 OE+OO,  1 . 00000000000E+00 ) , 


while  the  subroutine  (5.6)  for  CDIV  overflows  when  trying  to  compute  DEN DM. 

c.  Ease  of  use.  For  VAR  U,V,Wi  COMPLEX,  the  use  of  (5.6)  allows  one  to  write 


(5.13) 


V  V  /  W; 


in  the  source  code  for  the  program  instead  of 


(5.14) 


CDIV(V,W)» 


U  :■ 


as  in  ordinary  Pascal.  In  the  case  of  complicated  expressions  involving  complex  numbers, 
the  gain  in  programming  ease  using  ordinary  mathematical  notation  with  operators  instead  of 
function  and  procedure  calls  is  significant.  The  source  code  is  more  likely  to  be  correct 
in  the  first  place,  and  also  will  be  easier  to  document  and  read  later. 

d.  Compilation  time.  The  function  CD IV  has  to  be  compiled  from  the  source  code  (5.5) 
for  each  ordinary  Pascal  program  which  uses  complex  division.  The  operator  (5.6),  on  the 
other  hand,  is  given  by  pretranslated  code  which  is  automatically  linked  to  the  user's 
program  in  the  last  stage  of  the  compilation.  This  saves  a  considerable  amount  of 
compilation  time. 

The  accurate  complex  divison  in  Pascal-SC  turns  out  to  be  slower  than  the  function 
CDIV.  According  to  the  table  given  in  [4],  p.  267,  a  typical  time  for  the  complex  division 
U/V  is  100  milliseconds  for  a  2.5MHr  Z80  processor.  The  function  CDIV(U,V)  uses  six  real 
multiplications,  two  real  divisions,  and  three  real  addition/subtractions.  The  typical 
times  for  these  operations  given  in  (4]  total  60.4  milliseconds.  It  will  be  seen  later 
that  some  Pascal-SC  operations  are  actually  faster  than  their  inaccurate  real  simulations! 
however,  in  the  case  of  complex  multiplication  and  division,  one  pays  a  little  for 
guaranteed,  reliable  accuracy. 

In  addition  to  the  arithmetic  operators  +,  -,  *,  /  for  type  COMPLEX,  a  number  of 
additional  operators,  functions,  and  procedures  are  provided  in  the  Pascal-SC  complex 
library  for  convenience.  For  details  on  these,  including  the  domains  and  ranges  of  the 
standard  functions,  see  (28]. 

Sounded  complex  operations  +<,  +>,  -< ,  ->,  *<,  *>»  /<,  />  are  also  included  in  an 
external  library  in  the  form  of  pretranslated  code  [28].  Here,  rounding  is  carried  out 
componentwise .  Each  complex  number  s  “  (x,y)  with  |x|,|y|  <.  MAXREAL  will  belong  to  a 
rectangle  with  corners  which  are  the  floating-point  complex  numbers  A  »  (u,v),  B  » 

(u  ♦  4,  v),  C  «  (u  +  S,  v  +  n),  D  «  (u,  v  +  n),  and  which  contains  no  other  floating-point 


coaplex  numbers  (aee  Figure  5.1).  The  result  of  rounding  z  downward  will  be  V*  «  A  “ 

.(u,v),  while  z  Is  rounded  upward  to  It  *  C  ■  (u  +  t,  v  +  n),  where  4  and  n  are  the  spacinga 
In  the  floating-point  screen  In  the  horizontal  and  vertical  directions  In  the  coaplex 
plane,  respectively,  *nie  BPA  rounding  of  z  will  be  (BPA(x),  BPA(y)),  and  thus  could  be  any 
one  of  the  four  points  A,B,C,D.  For  example,  suppose 

(5.15)  z  -  (100  -  4i)/( 565  +  7891), 

which  is  not  a  coaplex  floating-point  number.  Paacal-SC  operations  give 

A  -  ( 100,-4 )/<( 565,789 )  -  ( 5.66437234668E-02,-8.6180350 1 157E-02 ) , 

(5.16) 

C  -  ( 100, -4)/> (565, 789 )  -  (5.66437234669E-02,-8.6180350ll56E-02) , 

while  the  BPA  for  (5.15)  la 

(5.17)  D  -  (100,-4 )/(565,789 )  -  (5. 66437234668E-02, -8.61803501 156E-02) . 


Figure  5.1.  A  Detail  of  the  COMPLEX  Floating-Point  Screen 


In  this  csss,  5  “  n  “  1*0  x  10~13,  and  the  rasults  (5.16)  and  (5.17)  locate  z  in  the 
rectangular  complex  interval  A  +  [0,  5.0*10  '*)  +  1(5.0x10”'*,  I.OxlO”'2),  according  to  the 
rules  f or  Fascsl-SC  rounding  (see  Figure  5.1). 

Complex  arithmetic  with  directed  rounding  can  be  used  as  the  basis  for  complex 
interval  arithmetic  (14). 

6.  Floating-point  INTERVAL  arithmetic.  Interval  arithmetic  (1)  [17],  (16)  is  based 
on  the  use  of  closed,  finite  intervals  (a,b)  of  real  numbers  as  its  basic  elements. 

Interval  arithmetic  has  a  number  of  significant  applications  in  scientific,  engineering, 
and  statistical  computation;  however,  its  use  has  not  been  widespread  up  to  now  because  of 
the  limitations  of  conventional  computer  arithmetic  units  (18)  and  ordinary  programming 
languages.  In  Pascal-SC,  interval  arithmetic  has  been  implemented  efficiently,  and  is  just 
as  convenient  to  use  as  real  or  complex  arithmetic. 

In  the  general  theory  of  computer  arithmetic  [14],  interval  arithmetic  is  regarded  as 
a  special  case  of  arithmetic  on  subsets  of  real  numbers.  Here,  if  X,  Y  are  subsets  of  R, 
and  o  e  {+,  -,  *,  /),  then 

(6.1)  Z  -  X  o  y  -  {xoy  |  x  €  X,  y  e  Y), 

by  definition.  For  division,  of  course,  0  €  Y  is  excluded.  If  X  -  [a,b]  and 
Y  “  [c,d]  are  intervals,  then  Z  ■  (r,s)  is  an  interval  if  defined,  and  the  endpoints  r,s  of 
Z  can  be  calculated  from  the  endpoints  of  X,Y  (14),  (17),  (16).  It  is  assumed,  of  course, 
that  the  intervals  X  ”  [a,b]  considered  are  finite  and  proper,  that  is,  a  <_  b. 

The  meet  basic  application  of  Interval  arithmetic  is  the  following;  If  the  value  of  a 
function,  for  example, 

(6.2)  w  -  9x4  -  y4  +  2y2, 


is  calculated  in  interval  arithmetic  for  intervals  X,Y,  then  the  resulting  interval  W  will 
contain  all  values  w  of  the  function  (6.2)  for  all  values  of  x  €  X,  y  €  Y.  This  property 


of  interval  arithmetic  allows  one  to  bound  the  ranges  of  functions  without  detailed 
analysis  of  maxima  and  minima.  Since  rounding  of  interval  operation*  is  outward  to  the 
smallest  floating-point  interval  which  contains  the  exact  reaulta  [28] ,  automatic  bounds 
for  round-off  error  can  be  obtained  conveniently,  by  taking  the  input  intervals  X,Y  to  be 
tingle  points,  that  is,  X  -  [x,x],  Y  «  ty,y].  Ml,  [17],  [16].  For  a  simple  application  of 
this  principle,  consider  the  evaluation  of  (6.2)  by  the  expressions 

(a)  w  s -  9*x»x*x*x  -  y*y*y*y  +  2*y*y» 

(6.3)  (b)  w  i-  9*(x**4)  -  y**4  +  2«(y**2)» 

(c)  w  t-  (3*(x**2)  -  y**2)»(3»(x*»2)  +  y**2)  +  2*(y**2>» 

in  REAL  arithmetic,  and 

(a)  W  s-  9*X*X*X*X  -  Y*Y*Y*Y  +  2*Y*Y» 

(6.4)  (b)  W  9*(X**4 )  -  Y**4  +  2*(Y**2 )i 

(c)  W  (3* (X**2 )  -  Y**2)*(3*(X**2)  +  Y**2)  +  2*(Y**2)» 

in  INTERVAL  arithmetic,  respectively,  the  power  operators  *‘  in  (6.3)(b)  and  (6.41(b)  are 
given  in  Appendix  A  in  the  source  code  for  the  program  WEVAL  which  was  used  to  calculate 
the  following  results. 

For  x  -  10864,  y  -  18817,  the  statements  (6.3)  give 


(a) 

W  “ 

1.58978  x  105, 

W  - 

[-1.841022  x  106, 

1.158978  x  106] 

(b) 

w  “ 

-8.41022  *  105, 

W  - 

(-8.410220  x  )05, 

1.158978  x  106] 

(c) 

w  * 

1.00000000000, 

V  - 

[  1.00000000000  , 

1.00000000000  ] 

The  enormous  width  of  M  in  (6.5)<a)  and  (6.5)(b)  indicates  that  the  values  given  for  w  can 
be  subject  to  large  roundoff  error,  and  hence  are  untrustworthy.  On  the  other  hand, 
evaluation  of  (6.2)  by  the  statements  (6.3)(c)  and  (6.4)(c),  respectively,  shows  that  this 


formulation  la  highly  accurate,  and  in  fact  proves  that  (6.3)(c)  yields  the  exact  value 
w  ■  1  of  the  function  (6.2)  for  the  given  values  of  x  and  y.  This  example  also  shows  that 
the  way  in  which  arithmetic  expressions  are  written  can  be  crucial  for  accuracy.  The 
techniques  for  evaluation  of  arithmetic  expressions  with  maximum  accuracy  [6]  can  be 
automated,  and  are  incorporated  in  a  Pascal-SC  demonstration  program  [26]. 

In  Pascal-SC,  floating-point  intervals  are  declared  in  the  stereotyped  way: 

TYPE  INTERVAL  -  RECORD  INF, SUP >  REAL  END I 

As  in  the  case  of  complex  numbers,  the  operators,  functions,  and  procedures  in  an  external 
library  manipulate  intervals  declared  in  this  form.  In  order  to  preserve  the  inclusion 
property  of  interval  arithmetic,  all  rounding  in  floating-point  INTERVAL  arithmetic  is 
outward:  The  result  given  for  xoy  is  the  smallest  floating-point  interval  Z  which  contains 
the  actual  result.  This  can  be  implemented  by  the  use  of  directed  rounding,  for  example, 
addition  can  be  performed  by  the  operator 

(6.7)  OPERATOR  +  ( A,B:  INTERVAL)  RES:  INTERVAL; 

VAR  C:  INTERVAL; 

BEGIN 

C.INF  «■  A. INF  ♦<  B.INF; 

C.SUP  s*>  A. SUP  +>  B.SUPf 
RES  I-  C; 

END; 

This  operator  is  available  as  a  pretranslated  subroutine,  declared  by  . 

OPERATOR  ♦  (A,Bl  INTERVAL)  RES:  INTERVAL 


EXTERNAL  68; 


If  I  denotes  a  generic  variable  of  type  INTERVAL,  and  K  of  type  INTEGER,  the  following 
15  arithmetic  operations  are  provided  for  interval  arithmetic: 

±1,  K  ±  I,  I  i  K,  I  ±  I, 

K  *  I,  I  *  R,  I  *  K, 

K  /  I,  I  /  R,  I  /  I. 


Because  REAL  floating-point  expressions  do  not  necessarily  yield  the  exact  values  of 
their  real  results,  operators  between  types  REAL  and  INTERVAL  are  not  included,  since  the 
resulting  interval  might  not  contain  the  true  outcome  cf  a  computation.  In  addition  to 
interval  versions  of  standard  functions  and  various  utility  procedures  [28],  the  interval 
library  includes  some  operators  which  work  with  intervals  as  sets  of  real  numbers.  There 
are  the  "lattice  operators" 


(6.8) 


Intersection, 
Interval  Hull, 


and  the  relational  operators 


<* 

Subinterval, 

>« 

Superinterval, 

>< 

Dis jointness, 

IN 

Point  Inclusion 

(28).  The  intersection  operator  **  will  generate  an  error  interrupt  if  its  operands  are 
disjoint  Intervals,  otherwise,  their  intersection  IOJ  will  be  computed.  If  I  *  [a,b]  and 
J  «  lc,d],  then  I  ++  J  ■  (min{a,c) ,max(b,d)]  is  the  smallest  interval  which  contains  both  I 
and  J.  The  relation  I  <-  J  is  TRUE  if  I  C  J,  otherwise  FALSE:  similarly,  I  >-  J  is  TRUE  if 


I  3  J.  The  result  of  I  X  J  is  TRUE  if  I  and  J  are  disjoint  intervals.  This  test  can  be 


used  to  avoid  I  **  J  in  this  case.  ..f  R  is  a  floating-point  number  (type  REAL),  R  IN  I  is 
TRUE  if  R  €  I  as  a  raal  number. 

With  regard  to  the  efficiency  of  the  implementation  of  interval  arithawtic  in  Pascal- 
SC,  Bohlender  and  Grflner  [4]  give  the  following  typical  times  in  milliseconds  for  a 
microcomputer  with  a  2.5MH*  Z80  processor i 

Operation  ♦  -  *  / 

REAL  2.2  2.2  6.0  10.0 

INTERVAL  5.4  5.4  23.0  31.0 


Considering  the  fact  that  each  interval  operation  has  to  determine  two  real  numbers,  and 
that  Interval  multiplication  and  division  require  the  calculation  of  four  real  products  in 
one  case  [14],  the  above  Indicates  an  almost  optimal  implementation  in  software.  By 
contrast,  factors  of  100  or  200  between  REAL  and  INTERVAL  arithmetic  speed  have  been  noted 
on  conventional  computers,  such  as  the  0K1VAC  1100  series.  [18].  In  addition,  while  a  REAL 
computation  provides  only  a  floating-point  number  which  approximates  the  result,  an 
INTERVAL  computation  on  the  other  hand  provides  an  Interval  which  is  guaranteed  to  contain 
the  true  result.  An  important  application  of  this  property  of  interval  calculations  will 
be  described  below  in  connection  with  the  solution  of  linear  systems  of  equations  with 
guaranteed  error  bounds. 

7.  Real  floating-point  vector  and  matrix  arithmetic.  Calculations  with  real  vectors 
and  matrices  are  among  the  most  commonly  encountered  tasks  in  scientific,  engineering,  and 
statistical  computation.  Pascal-8C  offers  the  user  the  same  reliability,  accuracy, 
controllability,  and  convenience  when  calculating  with  real  n-dimensional  floating-point 
vectors  and  matrices  as  it  does  for  the  scalar  types  REAL,  COMPLEX,  and  INTERVAL.  This  is 
accomplished  with  the  aid  of  an  external  source  code  library  KRLIB  of  operators,  functions, 
and  procedures  for  vector  and  matrix  manipulation,  and  a  built-in  function  SCALP  for  the 


calculation  of  acalar  products  of  f loating-point  vectors  to  BPA  accuracy  or  with  directed 
rounding  at  tha  option  of  tha  user* 

7.1.  Convenience.  In  order  to  illustrate  the  convenience  of  Pascal-SC  for  vector  and 
matrix  calculations,  suppose  that  A,  B  are  nxn  matrices,  and  x,  y,  z  are  n-dimensional 
vectors.  To  evaluate 

(7.1)  z  -  5.5ABX  3y , 

the  corresponding  expression  is  Paacel-SC  is 

(7.2)  z  i«  5. 5*A*B*X  +  3*y; 

which  uses  ordinary  operator  notation  instead  of  the  function  and  procedure  calls  which 
would  be  required  in  ordinary  Pascal  and  most  other  languages.  In  order  to  make  use  of  the 
software  provided  in  the  corresponding  external  library,  a  stereotyped  declaration  of 
floating-point  vector  and  matrix  data  types  is  expected* 

CONST  DIM  ■  I*  (The  actual  dimension  replaces  #} 

TYPE  DIMTYPE  -  1.  .DIM; 

R VECTOR  -  ARRAY  [DIMTYPE]  OP  REAL; 

RMATRIX  -  ARRAY  [DIMTYPE]  0 T  R VECTOR) 

The  operators  ♦,  -,  *,  and  the  operators  +<,  +>,  -<,  ->,  *<,  *>  with  directed  rounding 
are  available  for  various  permissible  combinations  of  opsrands,  for  example,  multiplication 
of  an  KVECTOR  by  an  INTEGER  or  REAL,  and  so  on  [28]. 

7.2.  Reliability,  accuracy,  and  controllability i  The  scalar  product  SCALP.  The 
general  theory  of  computer  arithmetic  [14]  requires  that  each  component  of  the  result  of  a 
vector  or  matrix  operation  be  rounded  to  the  BPA  for  the  actual  real  result,  or  downward  or 
upward  to  tha  closest  neighboring  floating-point  number  if  desired.  Addition  and 


subtraction  of  vectors  and  matrices  present  no  problems  from  the  standpoint  of  this 
requirement/  since  the  desired  results  can  be  calculated  componentwise  with  the  aid  of  the 
six  SEAL  arithmetic  operators  ±,  ±<,  ±>  described  in  {3.  Calculation  of  the  scalar 
products  of  vectors,  which  is  an  inherent  component  of  matrix  and  matrix  by  vector 
multiplication,  is  a  different  matter.  Ordinarily,  this  calculation  is  simulated  by  a  FOR 
loop  of  real  operations,  such  as  in  the  following  function: 

(7.3)  FUNCTION  SPROD(A,B:  R VECTOR ) !  REAL) 

VAR  I:  DIMTYPE;  S:  REAL; 

BEGIN 
S  Of 

FOR  1 1  —  1  TO  DIM  DO 
S  S  +  A(I)  *  B[I]» 

SPROD  S 
END) 

This  is  an  example  of  what  Kulisch  calls  the  •vertical*  definition  of  computer  arithmetic 
[14],  [12].  Of  course,  there  is  no  hope  that  the  result  of  SPROD  will  be  accurate  in 
general.  For  this  reason,  the  internal  calculations  in  a  function  of  this  kind  are  often 
done  in  higher  precision  than  the  external  calculation.  While  this  is  a  often  a  great  help 
in  some  cases,  it  still  does  not  solve  the  accuracy  problem.  On  the  other  hand,  the 
Pascal-8C  function 

(7.4)  SCALP(A,Bt  RVECTORf  ROUND)  INTEGER)  f 

which  is  a  built-in  feature  of  the  coaqpiler,  will  calculate  the  value  of  the  exact  scalar 
product  of  A  and  B  to  an  adjacent  floating-point  number,  with  rounding  downward,  upward,  or 
to  the  BPA  controlled  by  the  value  of  the  parameter  ROUND  [28] .  This  reliability, 
accuracy,  and  controllability  is  required  by  the  general  theory  of  casqxiter  vector  and 


matrix  arithmetic  [14].  It  can  be  achieved  by  special  algorithms  [14],  or  the  provision  of 
a  sufficiently  long  accumulator.  In  the  microcomputer  version  of  Pascal-SC,  this  "long 
accumulator*  is  implemented  in  software  [4],  but  the  same  thing  can  be  done  in  hardware 
[13],  and  can  be  expected  to  be  a  feature  of  future  advanced  mainframe  computers. 

In  order  to  allow  the  accumulation  of  several  scalar  products,  the  parameter  ROUND  can 
also  inhibit  the  clearing  of  the  long  accumulator  before  the  product  is  calculated  [28] . 

The  corresponding  values  are  given  in  the  following  table: 


Rounding 

Clear  Long  Accumulator 

Inhibit  Clearing 

BPA 

ROUND  -  0 

ROUND  -  4 

Downward 

ROUND  -  -1 

ROUND  -  3 

Upward 

ROUND  -  1 

ROUND  -  5 

If  the  long  accumulator  is  not  being  cleared,  other  arithmetic  operations  are  not  permitted 
between  successive  calls  of  SCALP  [28] . 

The  use  of  SCALP  makes  it  possible  to  calculate  the  results  of  matrix  and  matrix  by 
vector  multiplications  to  the  closest  floating-point  numbers,  or  rounded  to  the  closest 
larger  or  smaller  neighboring  values  if  desired.  This  reliability,  accuracy,  and 
controllability  distinguishes  Pascal-SC  vector  and  matrix  arithmetic  from  traditional 
packages. 

Some  of  the  important  properties  of  SCALP  are  illustrated  by  the  following  examples. 
a.  Accuracy.  For 

(7.5)  A  -  <1099,  10’",  -10”),  B  -  (1,  1, 

the  value  of  SPROD(A.B)  is  0,  of  course,  while  SCALP(A,B,0)  gives  the  correct  answer 
10  .  Persons  who  believe  that  multiple  precision  will  solve  all  accuracy  problems  should 

determine  how  much  precision  is  required  on  their  machine  for  SPROD(A,B)  to  duplicate  this 
result.  Once  satisfied,  they  can  then  try  (7.6)  below.  Although  these  examples  are 
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extreme ,  exanples  can  be  given  for  which  SPROD  is  highly  inaccurate  for  vectors  of  the 
types  one  can  expect  to  encounter  in  actual  problems. 

b.  Speed.  It  turns  out  that  the  accurate  scalar  product  function  SCALP  is  faster 
than  the  corresponding  FOR  loop  in  SPROD  if  DIM  >  1.  The  typical  times  given  in  [4]  in 
milliseconds  for  a  2.5MHz  Z80  processor  are: 

SCALP  8  +  5. 5 (DIM ) , 

FOR  LOOP  1  +  9.6(DIM) . 

Since  all  the  matrix  and  matrix  by  vector  multiplication  routines  are  based  on  SCALP,  they 
can  be  expected  to  execute  faster  than  their  inaccurate  simulations  in  REAL  arithmetic. 
Furthermore,  the  current  implementation  of  Pascal-SC  handles  access  to  elements  of  arrays 
very  efficiently  by  means  of  an  "array  descriptor"  (111,  from  which  addresses  of  elements 
can  be  calculated  quickly. 

c.  Deferred  overflow.  Under  ordinary  circumstances,  SCALP  will  not  indicate  an 
overflow  unless  the  actual  real  result  is  outside  the  range  of  representable  floating-point 
numbers.  For  example,  for 

(7.6)  A  -  <1099,  10'",  -10">,  B  -  HO99,  1,  1099), 

SCALP (A,B,0 )  will  compute  the  correct  result  10”",  while  SPROD(A,B)  will  overflow  for 

1-1. 

d.  Ease  of  use.  The  function  SCALP  is  just  as  easy  to  use  as  SPROD,  and  requires  no 
additional  source  code  in  the  program,  since  it  is  part  of  the  Pascal-SC  system. 
Furthermore,  SCALP  is  more  versatile,  since  its  arguments  can  be  arbitrary  one-dimensional 
arrays  of  floating-point  numbers  of  the  same  length,  of  which  RVECTOR  is  only  a  special 
case.  When  called,  SCALP  consults  the  array  descriptors  of  its  arguments  (11)  to  determine 
if  they  are  in  fact  of  equal  length,  and  then  calculates  their  scalar  product  if  this  is 
true.  Thus,  SCALP  can  be  used  to  calculate  scalar  products  of  vectors  of  various 


dimensions  in  the  same  program. 

8.  Complex  and  interval  floating-point  vector  and  matrix  arithmetic.  The  basic  ideas 
here  are  generally  the  same  as  for  real  sector  and  matrix  arithmetic:  Convenience,  based 
on  the  uae  of  operator  notation  for  expreasions,  and  accuracy.  The  subroutines  in  MCLIB 
and  MILIB,  respectively,  expect  the  type  declarations 

TYPE  CVECTOR  -  ARRAY  [DIMTYPE]  OF  COMPLEX) 

C MATRIX  -  ARRAY  [DIMTYPE]  OF  CVECTOR) 

and 

TYPE  IVECTOR  -  ARRAY  [DIMTYPE]  OF  INTERVAL) 

I  MATRIX  -  ARRAY  [DIMTYPE]  OF  IVECTOR) 

in  the  corresponding  cases. 

For  accurate  scalar  products,  the  respective  functions 

FUNCTION  CSCALP  (VAR  A,B:  CVECTOR)  AKDIM:  INTEGER):  COMPLEX) 

EXTERNAL  188) 

and 

FUNCTION  I  SCALP  (VAR  A,B:  IVECTOR)  AKDIM:  INTEGER):  INTERVAL) 

EXTERNAL  88) 

are  provided.  The  functions  form  the  basis  of  the  subroutines  for  accurate  matrix  and 
matrix  by  vector  multiplication.  The  products  are  computed  from  the  first  AKDIM  components 
of  each  vector  argument.  This  permits  the  flexibility  of  using  vectors  of  various 
dimensions  AKDIM  <_  DIM  in  the  same  program.  CSCALP  computes  the  BPA  for  the  scalar  product 
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of  complex  vectors;  directed  rounding  is  not  provided  for  the  complex  scalar  product  or 
complex  matrix  and  matrix  by  vector  multiplications.  ISCALP  computes  the  smallest  interval 
which  contains  the  exact  result,  as  in  the  case  of  other  interval  operations  [28). 

9.  Solution  of  linear  systems  of  equations  and  matrix  Inversion.  These  are  problems 
which  arise  time  after  time  in  scientific,  engineering,  and  statistical  computation,  ttie 
Pascal-SC  system  subroutines  for  these  purposes,  which  use  the  accurate  scalar  product  and 
interval  arithmetic,  yield  results  which  are  far  more  exact  than  can  be  obtained  by 
ordinary  floating-point  arithmetic.  Furthermore,  guaranteed  error  bounds  are  given  for 
results,  so  the  reliability  of  the  computation  is  immediately  determinable.  These 
subroutines  will  accept  either  real  or  interval  vectors  and  matrices. 

The  basic  procedure  in  the  system  library  LGLLIB  for  the  solution  of  linear  systems  of 
equations  is  LGLP  (an  acronym  for  the  German  words  for  "linear  equations  solution 
program").  This  procedure  is  declared  by 

PROCEDURE  LGLP (DIM, AKDIM:  INTEGER;  VAR  A:  RHATRIX,  VAR  Bs  R VECTOR; 

VAR  Y:  IVECTOR); 

EXTERNAL  524; 

[28) .  The  purpose  of  this  procedure  is  to  solve  the  linear  system 

(9.1)  Ax  »  B 

with  coefficient  matrix  A  and  right-hand  side  B.  Instead  of  a  floating-point  approximation 
to  the  solution  x,  LGLP  calculates  an  interval  vector  Y  which,  if  proper,  contains  the 
exact  solution  x  of  (9.1),  and  proves  that  the  floating-point  matrix  A  is  a  nonsingular 
real  matrix  (27).  This  allows  one  to  determine  not  only  an  approximate  value  for  x,  but 
also  guaranteed  error  bounds  for  it  [23 J.  The  parameter  AKDIM  in  the  formal  parameter  list 
allows  one  to  solve  systems  of  size  smaller  than  DIM  if  desired;  only  the  first  AKDIM  rows 
and  columns  of  A  and  components  of  B  are  involved  in  the  calculation. 


Failure  of  LGLP  to  return  a  proper  interval  vector  Y  indicates  that  A  is  singular  or 


extremely  ill-conditioned.  In  this  case,  the  components  of  Y  will  be  set  equal  to  the 
i^roper  interval  [+1,-1].  A  test  should  be  made  for  this  condition  immediately  on  return 
from  LGLF,  since  all  interval  subroutines  expect  proper  intervals  as  data  (28). 

Thus,  LGLP  either  gives  a  solution  with  guaranteed  accuracy  or  an  error  indication. 

In  practice,  LGLP  has  been  observed  to  succeed  for  well-known  examples  of  badly  conditioned 
matrices,  such  as  Hilbert  matrices  [26],  (27).  In  Appendix  B,  a  simple  program  to  solve 
linear  equations  of  order  up  to  20  is  given,  together  with  its  application  to  a  system  of 
five  equations  in  five  unknowns  in  which  the  coefficient  matrix  has  a  condition  number 
larger  than  4  x  10 1®.  Inspection  of  the  resulting  interval  vector  Y  shows  that  LGLP  was 
able  to  solve  this  system  to  an  accuracy  of  one  unit  in  the  twelfth  significant  digit.  The 
same  5x5  system  defeated  a  standard  linear  equation  solver  on  a  VAX  11/780. 

When  one  is  solving  several  systems  with  the  same  matrix  but  different  right  Bides, 
considerable  time  can  be  saved  if  the  approximate  LU-decomposition  of  A  and  other 
preliminary  calculations  are  only  done  once.  The  procedure  LGLPR  is  available  for  thiB 
purpose.  Its  declaration  is 

PROCEDURE  LGLPR  (DIM , AKDIM :  INTEGER;  VAR  A;  RMATRIX;  VAR  B:  HVECTOR > 

NRS:  BOOLEAN;  VAR  R;  RMATRIX;  VAR  MB:  IMATRIX; 

VAR  Y:  RVECTOR ) ; 

EXTERNAL  522; 


[28].  The  first  call  of  this  procedure  is  with  NRS  *■  FALSE.  Matrices  needed  to  process 
subsequent  right  sides  will  be  computed  and  stored  as  R  and  MB.  Subsequently,  LGLPR  is 
called  with  NRS  «  TRUE  for  each  new  right  side. 

Similar  procedures  LGLI  and  LGLIR  are  available  for  the  case  that  the  consonants  of 
the  coefficient  matrix  A  and  right  side  B  are  intervals,  that  is,  A  is  of  type  IMATRIX  and 
B  is  of  type  IVECTOR.  The  interval  vector  Y  in  this  case  bounds  all  solutions  of  real 
systems  with  coefficient  matrices  belonging  to  A  and  right  sides  belonging  to  B.  This  can 
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be  helpful  in  case  where  the  data  are  subject  to  uncertainty* 

For  matrix  inversion,  the  subroutine 

PROCEDURE  INVP  (DIM, AKDIM)  INTEGER)  VAR  A:  RMATRIX)  VAR  Cs  I MATRIX ) ) 

EXTERNAL  526) 

will,  if  successful,  compute  an  interval  matrix  C  which  contains  the  inverse  of  the  real 
(point)  matrix  A.  Singularity  or  extreme  ill-condition  of  A  is  reported  in  the  same  way  as 
for  LGLP,  while  successful  calculation  of  C  proves  that  a  is  nonsingular,  as  before. 
Finally, 

PRODCEDURE  INVI  (DIM, AKDIM)  INTEGER)  VAR  A,Ct  IMATRIX)) 

EXTERNAL  527) 

is  used  for  inversion  of  interval-valued  matrices.  If  C  is  computed  successfully  as  an 
IMATRIX  of  proper  intervals,  then  C  contains  the  inverses  of  all  real  matrices  contained  in 
the  interval  matrix  A.  Return  of  C  with  all  components  equal  to  the  improper  Interval 
1+1, -1)  indicates  that  A  contains  at  least  one  singular  or  very  badly  conditioned  real 
matrix. 

On  a  microcomputer  with  64Kb  of  storage,  LGLP  and  LGLPR  are  limited  to  about  DIM  »  25 
or  less,  LGLI,  LGLIR,  and  INVP  to  DIM  *  20,  and  INVI  to  DIM  «  15.  For  these  relatively 
small  systems,  the  time  required  for  execution  seems  to  be  reasonable.  As  in  the  case  of 
the  rest  of  the  Pascal-SC  system,  source  code  for  declarations  and  pretranslated  code  for 
the  above  procedures  can  be  found  in  an  external  library. 

10.  Eigenvalues  and  eigenvectors.  The  Pascal-SC  system  provides  the  standard 
subroutine 

PROCEDURE  EIGEN  (DIM, AKDIM)  INTEGER)  VAR  A:  RMATRIX)  LAMBDA)  REAL) 

VAR  X)  RVECTOR)  VAR  ILAMBDA:  INTERVAL)  VAR  Y)  I  VECTOR) ) 

EXTERNAL  534) 
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for  the  calculation  of  guaranteed  interval  bounds  for  real  eigenvalues  and  vectors  of  real 


matrices  A,  in  particular,  symmetric  matrices.  This  is  another  type  of  calculation  which 
occurs  often  in  engineering  and  other  scientific  computation.  In  addition  to  the  actual 
dimension  AKDIM  and  the  matrix  A,  EIGEN  expects  floating-point  approximations  LAMBDA  and  X 
to  the  eigenvalue  and  eigenvector  of  interest,  or  at  least  values  with  which  to  start  the 
calculation.  If  successful,  the  interval  value  ILAMBDA  and  Interval  vector  Y  returned 
Include  an  exact  real  eigenvalue  and  eigenvector  of  the  floating-point  matrix  A,  and 
furthermore  guarantee  that  the  included  eigenvalue  is  of  multiplicity  one.  Hence,  EIGEN 
will  not  succeed  for  multiple  eigenvalues  [si],  [2B] .  In  case  of  failure,  EIGEN  will 
return  Improper  intervals  for  ILAMBDA  and  the  components  of  Y. 

11.  The  accurate  sum  of  n  floating-point  numbers.  Statistical  calculations,  in 
particular,  often  require  the  computation  of  the  sum  of  n  floating-point  numbers  and 
perhaps  also  their  squares, 

n  n 

(11.1)  8  -  Z  a  ,  T  -  Z  a  . 

1-1  i-1 

In  Pascal-SC,  it  is  possible  to  compute  the  BPA  for  S  and  T,  or  round  the  result  upward  or 
downward  to  the  closest  floating  point  number  by  taking  the  a^  as  components  of  an  RVEcfOR 
A  and  using  SCALP,  For  E  -  (1,1,1,...,1),  one  has 

(11.2)  S  -  SC  ALP  ( A,  E,  ROUND);  and  T  -  SCALP(A,A ,  ROUND )  ; 

with  the  desired  best-possible  result.  However,  in  the  case  of  the  sum  S,  the  standard 
Pascal-SC  subroutine 

FUNCTION  SUM  (VAR  Ai  RVECTOR;  AKDIM ,  ROUND  I  INTEGER) :  REAL; 

EXTERNAL  460; 
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is  also  provided.  This  function  performs  the  addition  of  ths  first  AKDIM  elements  of  A  to 
the  BPA  or  result  of  directed  rounding  of  the  BPA.  SUM,  like  SCALP,  uses  the  long 
accumulator,  and  the  values  of  ROUND  have  the  same  significance  as  given  in  |7  for  SCALP. 

It  is  thus  possible  to  call  SUM  again  without  clearing  the  long  accumulator,  to  allow 
independent  accumulation  of  partial  sums  without  loss  of  accuracy.  However,  no  other 
arithmetic  operations  are  allowed  between  successive  calls  to  SUM  [28). 

When  computing  sums  of  interval  numbers,  one  can  use 

(11.3)  IS  ■-  ISCALPdA,  IE,DIM ) ; 

where  IE  has  components  all  equal  to  [1,1].  For  the  interval  sum  of  squares,  however,  it 
is  preferable  to  form  the  interval  vector  IQA  with  components  equal  to  I8QR(A[I] )  for  I  “ 
1..DIM,  and  then  compute 

(11.4)  IT  I-  ISCALPt IQA, IE, DIM) ; 

rather  than  ISCALP(IA,IA,DIM),  for  the  reason  given  in  {6  about  the  preferability  of 
ISQR(X)  to  X*X  for  intervals. 

12.  Programming  in  Pascal-8C.  The  only  new  techniques  in  Pascal-SC  for  a  Pascal 
programmer  to  acquire  are  the  definition  and  use  of  operators.  Except  for  these,  there  is 
no  difference  between  Pascal  and  Pascal-SC  programming.  Therefore,  the  discussion  here 
will  focus  on  the  operator  concept  [19].  The  definition  of  an  operator  subroutine  is 
headed  by 

OPERATOR  <name>  (<formal  parameter  llst>)  <result  name>:  <result  type>  r 

The  code  following  this  heading  is  the  same  as  for  a  function  having  the  same  purpose.  The 
result  must  be  assigned  to  <result  name>  in  its  entirety  before  leaving  the  subroutine. 

For  example,  if  <result  name>  -  RES  is  of  type  INTERVAL,  one  must  calculate  an  Interval  U 
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and  make  tha  assignment  RES  i“  Ui  before  leaving  tha  aubroutina,  rather  than  calculating 
RES. DIP  and  RES. SUP  aaparataly.  Tha  formal  parameter  liat  conaiata  of  ona  or  two 
idantifiara  and  thair  types.  Thua,  oparatora  ara  aithar  unary  or  binary.  Sinca  oparatora 
occur  in  axpraaaion  atringa  {'infix*  notation),  thair  arguments  hava  to  ba  of  axpraaaion 
typa.  That  ia,  VAR  A,  ate.,  ia  not  allowed  in  tha  formal  parameter  liat.  Tha  examplea  of 
operators  given  in  fS  and  $6  and  below  can  bo  used  aa  models. 

Thera  are  two  ways  to  name  an  operator  in  Pascal-SC: 

(i)  By  rodafining  ('overloading')  ona  of  tha  standard  Pascal-SC  operator  symbols  for 
a  new  data  typa  or  typaa. 

(ii)  By  use  of  an  arbitrary  name  selected  by  the  user  which  conforms  to  the  ordinary 
rules  for  Idantifiara  in  Pascal  (10].  In  this  caae  (see  below),  the  priority  of  tha 
operator  also  has  to  be  declared. 

These  two  Mthoda  will  be  discussed  separately. 

12.1.  Overloading  atandard  operator  symbols.  This  is  tha  most  common  method  used  in 
scientific  and  enginaaring  computing  to  name  Pascal-SC  operators,  since  ona  usually  wishes 
to  follow  the  ordinary  mathematical  notation  encountered  in  the  forsulas  being  used.  The 
atandard  operator  symbols  in  Pascal-SC  ara,  in  order  of  decreasing  priority! 

Unary  operators! 

MOT,  *  (unary),  -  (unary) 

Multiplicative  (binary)  operators: 

•  /  DIV  MOD  AMD  •*  *>  •<  />  /< 

Additive  (binary)  operators: 

+  -♦>♦<  ->  -<  ♦♦  OR 

Relational  (binary)  operators: 

»<><->-<>  IN  >< 

The  fundamental  distinction  between  a  unary  and  a  binary  OPERATOR  is  that  the  formal 


parameter  liat  for  the  operator  contains  exactly  one  parameter  in  the  first  case,  and 


exactly  two  in  the  second,  and  these  are  the  only  possibilities.  An  overloaded  operator 
will  have  the  ease  priority  as  its  symbol  in  the  table  above.  In  the  case  of  +  and  -,  the 
parasMter  list  of  the  operator  heading  will  specify  whether  they  are  unary  (highest 
priority)  or  binary. 

One  convenience  of  Pascal-SC  that  is  immediately  apparent  is  that  one  can  define  *•  to 
perform  exponentiation  on  whatever  numerical  types  are  appropriate  for  the  application  at 
hand.  However,  this  should  be  done  with  care.  Some  good  methods  are  given  in  [7).  Source 
code  for  a  simple,  'repeated  squaring*  [22]  implementation  of 

OPERATOR  **  (R:  REAL)  K:  INTEGER)  RES)  REAL; 

to  perform  R*  for  integral  powers  of  floating-point  numbers  is  given  in  Appendix  A.  This 
operator  makes  it  possible  to  write  x3,  x4,  etc.  as  x"3,  x"4,  etc.  in  expressions  to  be 
evaluated,  which  is  a  more  convenient  way  to  represent  these  simple  powers  than  by  a 
procedure  or  function  call  as  in  ordinary  Pascal.  (A  function  or  procedure  should  be  used 
to  compute  the  result  of  raising  an  Interval  base  to  an  interval  power,  since  the  operator 
**  is  used  to  compute  the  intersection  of  INTERVAL  variables  (see  (6.8)).) 

The  order  of  the  operands  in  the  formal  parameter  list  determines  the  order  in  which 
the  operator  will  be  applied.  The  compiler  distinguishes  various  uses  of  the  same  operator 
symbol  by  the  type(s)of  its  operand(s),  and  their  order  if  the  operator  is  binary.  Thus, 
in  the  same  program,  can  be  used  to  denote  addition  of  complex  numbers,  intervals, 
vectors,  matrices,  quaternions,  polynomials,  etc.,  in  addition  to  its  standard  meaning  for 
integers  and  floating-point  numbers.  All  that  is  required  is  that  the  appropriate 
definition  of  OPERATOR  +  be  given  in  the  heading  of  the  program  for  each  meaning  of  in 
the  body  of  the  program. 

Of  course,  the  compiler  recognizes  only  the  rules  of  arithmetic  for  user-defined  data 
types  which  are  provided  to  it  by  the  programmer.  For  example,  if  one  wishes  to  use 
expressions  in  which  variables  of  both  type  INTEGER  and  type  GRADIENT  [22]  appear,  both 
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OPERATOR  ♦  (Ki  INTEGER)  Gi  GRADIENT)  RESi  GRADIENT! 


•ltd 

OPERATOR  ♦  (Gi  GRADIENT!  Kt  INTEGER)  RE8 I  GRADIENT! 

•mat  be  defined  in  tha  haading  of  tha  program  ao  that  tha  compiler  can  produce  coda  for 
both  K  ♦  G  and  G  +  K.  Typa  GRADIENT  conaiata  of  tha  valua  of  a  function  together  with  ita 
gradient  vector,  and  ia  declared  by 

TYPE  GRADIENT  -  RECORD  Pi  READ!  DPi  RVECTOR  END| 

in  Pascal-SC  122).  In  thia  caaa,  both  oparatora  produce  the  sane  result,  consisting  of  tha 
alteration  of  the  function  value  G.P  of  tha  GRADIENT  variable  G  to  JC  +  G.p  -  g.P  *  X, 
respectively,  with  no  change  in  tha  gradient  vector  G.DP  [22].  However,  it  could  happen 
that  the  user  ia  working  with  quantities  for  which  addition  ia  not  necessarily 
coaaoutative.  Pascal-SC  allows  the  possibility  of  defining  tha  result  of  *+"  or  any  other 
binary  operator  to  ba  dependent  on  the  order  of  tha  operands. 

12.2.  Named  operators.  In  Pascal-SC,  the  user  can  name  operators  according  to  the 
ordinary  Pascal  rules  for  identifiers  (10).  For  example,  the  factorial  operator  (a  unary 
operator)  could  be  called  PAC.  In  this  case,  PAC  4  would  have  the  value  41  *  24.  Note 
that  parentheses  are  not  used  unless  the  operation  is  applied  to  an  expression.  Similarly, 
a  named  binary  operator  is  written  between  its  operands  ("infix'  notation)  in  the  same  way 
as  +,  -,  *,  /,  etc.  The  operator  FAC  can  be  defined  in  the  program  heading  as  followss 

OPERATOR  PAC  (Al  INTEGER)  RES  I  INTEGER! 

BEGIN 

IP  A  <-  1  THEN  RES  I*  1 

EDSE  RES  !•  A  •  PAC  (A  -  1)> 

ENDi  {  Recursive  definition  of  OPERATOR  FAC  } 
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In  terms  of  FAC,  the  binomial  coefficient  C(N,X)  could  be  computed  by  use  of  the  statement 

(12.1)  B1N0M  l-  FAC  N  DIV  (FAC  K  *  FAC  <N  -  K))| 

This  example  is  for  only  for  illustration  of  the  construction  of  a  named  operator  and  the 
possibility  of  recursion.  On  the  microcomputer  implementation  of  Pascal*-8C,  INTEGER 
arithmetic  is  implemented  only  for  integers  I  such  that  -32768  _<  I  _<  32767  [28],  and  thus 
FAC  N  can  be  computed  only  for  H  <  1.  Actual  computation  of  factorials  should  be 
accomplished  by  a  type  conversion  to  REAL,  and  controlled  by  WHILE  or  UNTIL,  in  order  to 
avoid  stacking  recursions  too  deeply. 

Similarly,  the  binary  Boolean  operator  XOR  for  "exclusive  or"  could  be  defined  by 

OPERATOR  XOR  (A.Bi  BOOLEAN)  RESi  BOOLEAN) 

BEGIN 

RES  (A  AND  NOT  B)  OR  (NOT  A  AND  B) I 
END) 

A  typical  program  statement  using  XOR  would  be 

(12.2)  IF  0BS1  XOR  OBS2  THEN  PROB  I"  0.25  ELSE  PROB  ««  0.75) 

which  would  assign  the  value  0.25  to  PROB  if  just  one  of  OBSl,OBS2  is  TRUE,  or  0.75 
otherwise. 

In  order  for  the  Pascal-SC  compiler  to  recognize  FAC  and  XOR  as  the  names  of 
operators,  and  assign  priorities  to  them,  a  PRIORITY  declaration  for  each  named  operator 
must  follow  directly  after  the  heading  line  of  the  program,  which  gives  the  name  of  the 
program  and  the  list  of  files  used.  From  highest  to  lowest  priority,  these  priority 
declarations  have  the  forms 


PRIORITY  <Op*rator  nui>  “  •>  {  (Mary  operator*  } 


PRIORITY  <Op«rator  name>  ■  *>  {  Multiplicative  operators  } 

PRIORITY  <Operator  na**>  •  +i  {  Additive  operators  > 

PRIORITY  <Operator  name>  “  ”i  {  Relational  operators  } 

Thus,  suppose  one  writes  a  program  called  CHANCE  which  uses  the  operators  FAC  and  XOR 
to  calculate  probabilities  of  outcomes  in  some  stochastic  model ,  and  only  the  standard 
files  INPUT  and  OUTPUT  are  used  for  communication  between  the  program  and  the  outside 
world.  If  XOR  is  to  have  the  same  priority  as  OR,  then  the  first  three  lines  in  the 
heading  of  the  source  code  for  the  program  would  be 

PROGRAM  CHANCE  (INPUT, OUTPUT); 

PRIORITY  FAC  -  Cl 
XOR  -  +| 

The  standard  sequence  of  definitions  and  declarations  would  then  follow  to  complete 
the  heading  of  the  program,  and  then  the  body  of  the  program  consisting  of  the  actual 
statements  to  be  executed.  The  structure  of  a  Pascal-SC  program  therefore  differs  only 
slightly  from  that  of  an  ordinary  Pascal  program,  as  shown  in  the  next  section. 

12,3.  Structure  of  a  Paacal-SC  program.  In  Pascal-SC,  the  order  of  declarations  in 
the  heading  la  somewhat  freer  than  in  standard  Pascal  [19].  However,  the  general  principle 
applies  that  everything  must  be  declared  or  defined  before  use.  Since  overloading  standard 
operator  symbols  is  more  common  than  using  named  operators,  the  headings  of  most  Pascal-SC 
programs  will  look  identical  to  Pascal  programs  except  for  the  OPERATOR  definitions. 
Programming  is  further  simplified  because  Pascal-SC  already  provides  operators  for  most  of 
numerical  data  types  commonly  encountered  in  scientific  and  engineering  computation,  such 
as  complex  numbers,  intervals,  vectors,  and  matrices,  as  explained  in  the  previous 
sections. 

The  difference  between  Pascal  and  Pascal-SC  programs  is  most  striking  in  the 


-32- 


statements  of  the  actual  body  of  the  program.  Here,  the  power  of  operator  notation  makes 
la  poaaible  to  write  expressions  clearly  and  compactly.  This  elimination  of  complicated 
sequences  of  function  and  procedure  calls  shortens  programs  and  makes  the  source  code  much 
easier  to  read  and  understand.  This  facilitates  documentation  as  well  as  use  of  the 
program. 

With  just  two  exceptions,  noted  below  in  boldface  type,  the  sequence  of  a  Pascal-SC 
program  is  identical  to  an  ordinary  Pascal  program  [10]: 

PROGRAM  <Mame>  (  CList  of  internal  file  names>  )< 

PRIORI  TT 
LABEL 
CONST 
TYPE 

VAR  declarations  and  definitions  of  the  program  heading> 

PROCEDURE 

FUNCTION 

OPERATOR 

begin 

. .  <Statements  comprising  the  body  of  the  program> 

END. 

The  order  in  which  procedures,  functions,  and  operators  are  declared  is  arbitrary,  as 
in  Pascal.  The  implementation  of  the  operator  concept  in  Pascal-SC  is  based  on  the  fact 
that  the  underlying  virtual  machine  (the  so-called  KL/F  machine)  can  stack  operands  of 
arbitray  data  types,  so  that  functions  with  results  of  arbitrary  type  can  be  computed 
efficiently  [11].  This  refinement  is  permits  considerable  savings  in  the  number  of  machine 
instructions  actually  needed,  and  hence  leads  to  shorter  execution  times  [11]. 

13.  Conclusions.  The  accuracy  of  Pascal-SC  arithmetic  and  the  convenience  of 
operator  notation  for  manipulation  of  numerical  and  other  data  types  make  this  language  a 
valuable  tool  for  scientific,  engineering,  and  statistical  computation.  In  addition  to  its 
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usefulness  for  routine  problems,  such  ss  solution  of  linear  systems  of  equations, 
experience  has  shown  that  it  is  possible  use  this  language  to  program  and  carry  out  some 
rather  sophisticated  computations,  even  on  a  microcomputer.  Examples  include  numerical 
solution  of  a  nonlinear  integral  equation  [25],  the  solution  of  nonlinear  systems  of 
equations  by  iterative  methods  [22],  [9],  and  the  solution  of  ordinary  differential 
equations  by  real  and  interval  Taylor  series  [8].  In  these  applications  the  operator 
concept  of  Pascal-SC  was  used  to  implement  automatic  evaluation  of  derivatives  and  Taylor 
series  for  functions  defined  by  expressions  in  ordinary  mathematical  notation  [8],  [22], 
[24].  The  microcomputer  systems  used  in  these  investigations  can  best  be  described  as 
minimal:  Eight-bit  machines  with  Z80  processors,  64Kb  of  main  storage,  and  two  disk 
drives.  The  Pascal-SC  compiler  used  was  developed  by  Profs.  U.  Kulisch  and  and  H.-W. 
wlppermann  and  their  associates  at  the  Universities  of  Karlsruhe  and  Kaiserslautern  in 
Germany,  and  is  described  in  [3],  [19],  and  [28].  Even  these  modest  resources  appear 
adequate  for  many  of  the  day-to-day  calculations  needed  by  engineers,  scientists,  and 
statisticians,  as  well  as  for  research  on  methods  in  numerical  analysis  which  can  be 
applied  to  larger  problems.  As  "personal  computers"  grow  in  size  and  speed,  the  accuracy 
and  convenience  of  Pascal-SC  will  provide  the  user  with  a  more  powerful  tool,  and  its 
features  will  also  be  advantageous  on  forthcoming  larger  machines. 

14,  Acknowledgments.  The  author  would  like  to  thank  Professor  George  F.  Corliss  for 
reading  the  manuscript  carefully,  and  a  large  number  of  valuable  suggestions.  Example 
(6.2)  was  provided  by  Prof.  Dr.  Ulrich  Kulisch,  and  the  system  of  equations  solved  by  the 
Pascal-SC  program  in  Appendix  B  was  one  of  several  similar  systems  solved  for  Dr.  S.  Takagi 
of  the  U.  S.  Army  Cold  Regions  Research  and  Engineering  Laboratory. 
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Evaluation  of  (6.2)  w  =  9x 
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in  Pascal-SC 


1.  Source  code  for  the  program  WEVAL. 


PROGRAM  VEVAL( INPUT, OUTPUT); 

(*  This  program  calculates  v  “  3*(x**4)  -  y**4  +  2*(y**2)  in  real  and 
interval  arithmetic.  *) 

|USES  INTERVAL;  (*  DIRECTS  COMPILER  TO  USE  INTERVAL  LIBRARY  *) 

VAR  C:  CHAR;  w,x,y:  REAL;  V,X,r:  INTERVAL; 

(*  POWER  OPERATORS  *) 

OPERATOR  **  (Rs  REAL ; K :  INTEGER)  RES:  REAL;  (*  R  **  K  *) 

VAR  L:  INTEGER jU:  REAL; 

BEGIN  (*  OPERATOR  R  **  K  *) 

IF  (R  -  0)  AND  (K  <-  0)  THEN 
BEGIN  (•  ERROR  *) 

WRITELNC EXPONENTIATION  ERROR,  0  **  K,  K  <-  O'); 

SVR(O)  (*  RETURN  TO  OPERATING  SYSTEM  *) 

END;  (*  ERROR  *) 

IF  (K  -  0)  OR  (R  -  1)  THEN  U:-l 

ELSE  IF  K  -  1  THEN  U:-R 

ELSE  (*  K  <>  0,1  *) 

BEGIN  (*  REPEATED  SQUARING  *) 

L:-ABS(K);U:-1; 

REPEAT 

IF  L  MOD  2  -  1  THEN  Ui-R*U; 

L:-L  DIV  2; 

IF  L  <>  0  THEN  R:-R*R 
UNTIL  L  -  0; 

IF  K  <  0  THEN  U:-l/U  (*  NEGATIVE  EXPONENT  *) 

END;  (*  REPEATED  SQUARING  *) 

RESi-U 

END;  (*  OPERATOR  R  **  K  *) 


OPERATOR  **  (I i  INTERVAL;*:  INTEGER)  RES:  INTERVAL;  <*!****) 


VAR  L:  INTEGER; U:  INTERVAL; 

BEGIN  (*  OPERATOR  I  **  K  *) 

IP  (0  IN  I)  AND  (K  <-  0)  THEN 
BEGIN  (*  ERROR  *) 

WRITELNC EXPONENTIATION  ERROR,  I  **  K,  0  IB  I’)l 
SVR(O)  (*  RETURN  TO  OPERATING  SYSTEM  *) 

END;  (*  ERROR  *} 

IF  K  -  0  THEN 

BEGIN  (*  K  -  0  *) 

U.INF:-1;U.SUP:-1 
END  (*  K  -  0  *) 

ELSE  IF  K  «  1  THEN  U.— I 

ELSE  (*  K  <>  0,1  *) 

BEGIN  (A  REPEATED  SQUARING  *) 

L:-ABS(K);U.INF:-1;U.SUP:-1; 

REPEAT 

IF  L  MOD  2  -  1  THEN  U:-I*U; 

L:-L  DIV  2; 

IF  L  <>  0  THEN  I:«ISQR(I); 

UNTIL  L  -  0; 

IF  K  <  0  THEN  U:-l/U  (*  NEGATIVE  EXPONENT  *) 

END;  (*  REPEATED  SQUARING  *) 

RES:-U 

END;  (*  OPERATOR  R  **  K  *) 

(*  END  OF  POWER  OPERATORS  *) 

PROCEDURE  WWRITE(w:  REAL;  W:  INTERVAL); 

BEGIN  (*  WRITE  RESULTS  OF  REAL  AND  INTERVAL  EXPRESSIONS  *) 

WRITELN;WRITELN( '  w  -  ',»);  (*  OUTPUT  OF  RESULTS  *) 

WRITELN;WRITELN( '  W-  [ ' .W.INF , ’ , ' ,W.SUP,' j ' ) ; 

END;  (*  WRITE  RESULTS  *) 


BEGIN  (*  MAIN  PROGRAM  *) 


C  s-  *Y' ;  WHILE  C  -  'Y'  DO 
BEGIN  (*  ACTUAL  CALCULATION  *) 

WRITELNC 'ENTER  VALUES  OP  X  AND  Y'  )  jREAD(x,y)  ; 

X  s-  INTPTCx) ;  Y  t-  INTPTCy);  (*  CONVERT  REAL  VALUES  TO  INTERVALS  *) 

»  : -  9*x*x*x*x  -  y*y*y*y  +  2*y*y;  (*  REAL  RESULT  *) 

V  i-  9*X*X*X*X  -  Y*Y*Y*Y  +  2*Y*Y;  (*  INTERVAL  RESULT  *) 

WRITELN ; 

WRITELN( ' (*)  w  «-  9*x*x*x*x  -  y+y+y+y  ♦  2*y*y;  glvaas'); 

UURITE(w ,U) j  (*  OUTPUT  RESULTS  *) 

w  s-  9*Cx**4)  -  y**4  +  2*(y**2)j  (*  REAL  RESULT  *) 

W  f  9*(X**4)  -  Y**4  +  2*(Y**2);  (*  INTERVAL  RESULT  *) 

WRITELN; 

URITELNC(b)  w  «-  9*(x**4)  -  y**4  ♦  2*(y**2>;  giva.i'); 

WWRITE(«,W);  (A  OUTPUT  RESULTS  *> 

»  (3*(x**2)  -  y**2)*(3*(x**2)  «•  y**2)  +  2*(y^2);  (*  REAL  RESULT  *) 

W  (3*(X**2)  -  Y**2)*(3*(X**2)  +  Y**Z)  +  2*(Y**2); 

(*  INTERVAL  RESULT  *> 


WRITELN; 

WRITELNC ' (c)  w  (3*<x**2)  -  y**2)*C3*Cx**2)  +  y**l)  +  2*(y**2);  glvaai1); 
WWRITE(w.W);  (*  OUTPUT  RESULTS  *) 

WRITELN; 

WRITELNC 'MORE  VALUES  (Y/N)7 ' ) ;READ(C,C)  (*  CONTINUE  OR  QUIT  *) 

END  (*  ACTUAL  CALCULATION  *) 

END.  C*  MAIN  PROGRAM  *) 
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2.  Sample  results  using  the  program  WEVAL. 

B>XQPC  WEVAL 

ENTER  VALUES  OP  X  AND  Y 

*10864  18817 

(a)  w  i”  9*x*x*x*x  -  y*y*y*y  +  2*y*yj  gives: 

«  -  1 . 58978000000E+05 

W  -  ( -1 . 841 02  200000E+06 ,  1 . 15897800000E+06] 

(b)  v  i-  9*(x**4)  -  y**4  +  2*(y**2) j  gives: 
v  -  -  8 . 410 2  2000000E+0 5 

W  -  [ -8.41022000000E+05 ,  1.15897800000E+06] 

(c)  v  i-  (3*(x**2)  -  y**2)*(3*(x**2)  *  y**2)  +  2*(y**2);  gives: 

v  «  1 . OOOOOOOOOOOE+OO 

W  -  [  1. OOOOOOOOOOOE+OO,  1. OOOOOOOOOOOE+OO] 

MORE  VALUES  (Y/N)7 
*N 

KL/P-STOP 
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1.  Source  code  for  the  program  LINSYS. 


PROCRAM  LINSYS ( INPUT , OUTPUT ) j 

(*  Thla  program  coaputea  an  interval  vector  Y  containing  the  aolutlon 
of  linear  ayateaa  of  aquatlona  AX  -  B  up  to  order  20 »  or  returne  an 
arror  aaaaaga  if  the  coefficient  aatrlx  A  la  alngular  or  axtreaely 
badly  conditioned.  The  aatrix  A  ia  read  by  rowa,  followed  by  B.  *) 

loses  LCL, DIM-20;  (*  SETS  DIMENSION  FOR  EXTERNAL  LIBRARY  ROUTINES  A) 

VAR  A i  RKATRZXj  (*  COEFFICIENT  MATRIX  *> 

Bt  R VECTOR;  (*  RIGHT  SIDE  *) 

Yi  IVECTOR;  (A  INCLUSION  OF  SOLUTION  OF  AX  -  B  *) 

AKDIMi  INTEGER;  (*  ACTUAL  SIZE  OF  SYSTEM  A) 

I.Jl  DIMTYPE;  (A  INDEX  VARIABLES  *) 

C:  CHAR;  (A  CONTROL  VARUBLE  *) 

BEGIN  (A  MAIN  PROGRAM  *) 

Cl-'Y'i  WHILE  C  -  'Y'  DO 

BEGIN  (*  SOLUTION  OF  SYSTEM  *> 

WRITELN( ' INPUT  DIMENSION' ) ;READ( AKDIH) ;WRITELN; 

WRITELHC INPUT  MATRIX  BY  ROWS'); 

FOR  Ii«l  TO  AKDIM  DO  FOR  Jt-1  TO  AKDIH  DO  READ(A[I,J]) ; 

WRITELN; WRITELN( 'INPUT  RIGHT  SIDE'); 

FOR  Ii-l  TO  AKDIM  DO  READ(B(I}); 

LGLP(D1M, AKDIM, A,B,Y);  (*  SOLVE  SYSTEM  *) 

IF  Y [ 1 } . INF  <-  Y[ 1 J .SUP  THEN  (*  Y  IS  PROPER  *) 

BEGIN  (*  OUTPUT  OF  SYSTEM  AND  RESULTS  *> 

WRITELN; 

WRITELN( 'SOLVE  AY  -  B  WITH  INTERVAL  INCLUSION  OF  ANSWER'); 
WRITELN; 

FOR  Js-l  TO  AKDIM  DO  (*  OUTPUT  OF  A  BY  COLUMNS  *) 

BEGIN 

FOR  I i-l  TO  AKDIM  DO 

WRITELN('Al',I»2,',',J»2,’l  -  ' ,A[I,J]) (WRITELN; 

END;  (★  OUTPUT  OF  A  A) 

WRITELN; 

FOR  It-1  TO  AKDIM  DO 

WRITELN('B(',Ii2,'l  -  * (A  OUTPUT  B  *) 

WRITELN ;FOR  It-1  TO  AKDIM  DO  (*  OUTPUT  Y  A) 
WRITELN('Yl',Ii2,']  -  [ ' ,Y 1 1 ] .INF , ' , ' ,Y[I]  .SUP,  * J * ) ; 

END  (a  OUTPUT  OF  SYSTEM  AND  RESULTS  a) 
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ELSE  (*  Y  IS  IMPROPER  *) 

BEGIN  (*  ERROR  MESSAGE  *) 

VRITELN) VRITELN( 'THE  MATRIX  IS  SINGULAR  OR  BADLY  CONDITIONED’); 
FOR  t«-l  TO  AKDIM  DO  Y(I] .SUPl-Y[l] .INF  (*  RESET  Y  *) 

END;  (*  ERROR  MESSACE  *) 

VRITELN ;WRITELN( 'ENTER  ANOTHER  SYSTEM  (Y/N)?'); 

READ(C.C); 

END  (*  SOLUTION  OF  SYSTEM  *) 

END.  (*  MAIN  PROGRAM  *) 
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2.  Sample  results  using  LINSYS. 


A>XQP  LINS VS  TAX. DAT  CON» 

INPUT  DIMENSION 

INPUT  MATRIX  BY  ROWS 

INPUT  RIGHT  SIDE 

SOLVE  AY  -  B  WITH  INTERVAL  INCLUSION  OF  ANSWER 

A[  1,  1]  -  - 1. 91 56942 1219E+11 
A(  2,  1]  -  O.OOOOOOOOOOOE+OO 
A[  3,  1]  -  8.54599767833B+08 
A[  4,  1]  -  O.OOOOOOOOOOOE+OO 
A(  S,  lj  -  O.OOOOOOOOOOOE+OO 

At  1,  2]  -  O.OOOOOOOOOOOE+OO 
A(  2,  2]  -  6.12580328713E+U 
At  3,  2]  -  O.OOOOOOOOOOOE+OO 
A[  4,  2]  -  4.67810493936B+07 
At  5,  2]  -  2 . 73290749688E+09 

A(  1.  3]  -  l.OOOOOOOOOOOE+OO 
At  2,  3]  -  -l.OOOOOOOOOOOE+OO 
At  3,  3]  -  1.12080226343E+02 
At  4,  31  -  O.OOOOOOOOOOOE+OO 
A{  5,  3]  -  1.12075040461R+02 

A(  1,  4}  -  O.OOOOOOOOOOOE+OO 
A(  2,  4]  -  5. 21250000000E-01 
A(  3,  4]  -  O.OOOOOOOOOOOE+OO 
At  4,  4]  -  l.OOOOOOOOOOOE+OO 
At  5,  4]  •  O.OOOOOOOOOOOE+OO 

At  1,  5]  -  -5.03360090692E-03 
A(  2,  5]  -  1.60975964126E-02 
A(  3,  3]  -  2 . 2433360488 7E-03 
At  4,  5]  -  O.OOOOOOOOOOOE+OO 
A[  3,  3)  -  7 . 18104363616E-05 


B(  1]  -  1.83094646224E+04 
B[  2J  -  -5.91833428873E+04 
B(  31  -  -8.25723911920E+01 
B(  4]  -  1.9236966 7381E+33 
B(  3]  -  2.13844970037E+35 

Tt  X I  -  l  7.93624477937E+29,  7.93624477938E+29] 
Yt  21  -  (  7. 937390343998+29,  7.93739034600E+291 
Y[  3}  -  t '6.200037 60703E+32 , -4 . 200037607028+32 ] 
Yt  4]  -  ( -3. 71300190864E+37,-3.71300190863E+37) 
Y{  3]  -  (-3.020386104018+43,-3.020386104008+43) 

ENTER  ANOTHER  SYSTEM  (Y/N)t 
KL/P-STOP 
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3.  Contents  of  the  data  file  TAX. DAT. 


A>TYPE  TAK.DAT 
3 

- 1. 9136942121 9E+11  010  -5.03360090692E-03 
0  6.12380328713E+11  -1  5.2125E-01  1.60975964126E-02 
8.54399767833E+08  0  1.12080226545E+02  0  2 .2453S604687E-05 
0  4 . 67810493936E+07  010 

0  2 . 73290749688E+09  1.12075040461E+02  0  7 .18104363616E-05 
1 . 83094646224E+04  -5.91833428873E+04  -8.25723911920E+01  1.92589687381E+33 
2. 1 58449 70057 E+3 5 
N 
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ABSTRACT  (continued) 


(1)  Accurate  floating-point  arithmetic  for  real,  complex,  and  interval 
numbers,  vectors,  and  matrices.  Expressions  involving  these  types  can  be 
written  using  essentially  ordinary  mathematical  notation,  making  programming  and 
documentation  easier.  Numerical  results  are  generally  obtained  with  more 
accuracy  than  with  conventional  floating-point  arithmetic}  in  particular,  scalar 
products  of  vectors  are  computed  to  the  closest  float in g-point  number.  Standard 
routines  are  provided  which  return  guaranteed  error  bounds  as  well  as  answers 
for  the  solution  of  linear  systems  of  equations,  matrix  inversion,  and 
eigenvalue-vector  calculations. 

(ii)  User— defined  operators  to  allow  the  manipulation  of  data  of  various 
nonstandard  types  by  expressions  written  in  essentially  ordinary  mathematical 
notation,  instead  of  the  sequences  of  calls  to  functions  and  procedures  required 
in  ordinary  Pascal.  This  capability  permits  simple  and  straightforward  use  in 
programs  of  various  coordinate  systems  for  representation  of  variables, 
arithmetic  for  polynomials  and  series,  automatic  differentiation,  etc.,  and  thus 
adds  considerable  flexibility  and  power  to  the  language. 

Brief  descriptions  of  these  useful  features  of  PaBcal-SC  are  given, 
together  with  illustrative  examples.  Some  familiarity  with  ordinary  Pascal  is 
assumed. 


